Как дифференцировать матрицы и дифференцировать по матрицам: всё, что вам не рассказали про дифференцирование на матанализе
Любая задача машинного обучения — это задача оптимизации, а задачи оптимизации удобнее всего решать градиентными методами (если это возможно, конечно). Поэтому важно уметь находить производные всего, что попадается под руку. Казалось бы, в чём проблема: ведь дифференцирование — простая и понятная штука (чего не скажешь, например, об интегрировании). Зачем же как-то специально учиться дифференцировать матрицы?
Да в принципе-то никаких проблем: в этом параграфе вы не узнаете никаких секретных приёмов или впечатляющих теорем. Но, согласитесь, если исходная функция от вектора x имела вид f(x)=∣∣Ax−b∣∣2 (где A — константная матрица, а b — постоянный вектор), то хотелось бы уметь и производную выражать красиво и цельно через буквы A, x и b, не привлекая отдельные координаты Aij, xk и bs. Это не только эстетически приятно, но и благотворно сказывается на производительности наших вычислений: ведь матричные операции обычно очень эффективно оптимизированы в библиотеках, чего не скажешь о самописных циклах по i,j,k,s. И всё, что будет происходить дальше, преследует очень простую цель: научиться вычислять производные в удобном, векторно-матричном виде. А чтобы сделать это и не сойти с ума, мы должны ввести ясную систему обозначений, составляющую ядро техники матричного дифференцирования.
Основные обозначения
Вспомним определение производной для функции f:Rm→Rn. Функция f(x) дифференцируема в точке x0, если
f(x0+h)=f(x0)+[Dx0f](h)+oˉˉ(∣∣h∣∣),
где [Dx0f] — дифференциал функции f: линейное отображение из мира x-ов в мир значений f. Грубо говоря, он превращает «малое приращение h=Δx» в «малое приращение Δf» («малые» в том смысле, что на о-малое можно плюнуть):
f(x0+h)−f(x0)≈[Dx0f](h)
Отметим, что дифференциал зависит от точки x0, в которой он берётся: [Dx0f](h). Под ∣∣h∣∣ подразумевается норма вектора h, например корень из суммы квадратов координат (обычная евклидова длина).
Давайте рассмотрим несколько примеров и заодно разберёмся, какой вид может принимать выражение [Dx0f](h) в зависимости от формы x. Начнём со случаев, когда f — скалярная функция.
Примеры конкретных форм [Dx0f](h), когда f — скалярная функция
f(x) — скалярная функция, x — скаляр. Тогда
f(x0+h)−f(x0)≈f′(x0)⋅h
[Dx0f](h)=f′(x0)h=h⋅f′(x0)
Здесь h и f′(x0) — просто числа. В данном случае [Dx0f] — это обычная линейная функция.
f(x) — скалярная функция, x — вектор. Тогда
f(x0+h)−f(x0)≈i∑∂xi∂fx=x0hi,
то есть
[Dx0f](h)=(∇x0f)Th=⟨∇x0f,h⟩,
где ⟨∙,∙⟩ — операция скалярного произведения, а ∇x0f=(∂x1∂f,…,∂xn∂f) — градиент функции f.
f(x) — скалярная функция, X — матрица. Дифференциал скалярной функции по матричному аргументу определяется следующим образом:
f(X0+H)−f(X0)≈i,j∑∂Xij∂fX=X0Hij
Можно заметить, что это стандартное определение дифференциала функции многих переменных для случая, когда переменные — элементы матрицы X. Заметим также один интересный факт:
ij∑AijBij=trATB,
где A и B — произвольные матрицы одинакового размера. Объединяя оба приведённых выше факта, получаем:
Можно заметить, что здесь, по аналогии с примерами, где x — скаляр и где x — вектор (и f(x) — скалярная функция), получилось на самом деле скалярное произведение градиента функции f по переменным Xij и приращения. Этот градиент мы записали для удобства в виде матрицы с теми же размерами, что матрица X.
В примерах выше нам дважды пришлось столкнуться с давним знакомцем из матанализа: градиентом скалярной функции (у нескалярных функций градиента не бывает). Напомним, что градиент ∇x0f функции в точке x0 состоит из частных производных этой функции по всем координатам аргумента. При этом его обычно упаковывают в ту же форму, что и сам аргумент: если x — вектор-строка, то и градиент записывается вектор-строкой, а если x — матрица, то и градиент тоже будет матрицей того же размера. Это важно, потому что для осуществления градиентного спуска мы должны уметь прибавлять градиент к точке, в которой он посчитан.
Как мы уже имели возможность убедиться, для градиента скалярной функции f выполнено равенство
[Dx0f](x−x0)=⟨∇x0f,x−x0⟩,
где скалярное произведение — это сумма попарных произведений соответствующих координат (да-да, самое обыкновенное).
Посмотрим теперь, как выглядит дифференцирование для функций, которые на выходе выдают не скаляр, а что-то более сложное.
Примеры [Dx0f](h), где f — это вектор или матрица
Матрица, выписанная в предпоследней выкладке, — это знакомая вам из курса матанализа матрица Якоби.
Простые примеры и свойства матричного дифференцирования
Производная константы. Пусть f(x)=a. Тогда
f(x0+h)−f(x0)=0,
то есть [Dx0f] — это нулевое отображение. А если f — скалярная функция, то и ∇x0f=0.
Производная линейного отображения. Пусть f(x) — линейное отображение. Тогда
f(x0+h)−f(x0)=f(x0)+f(h)−f(x0)=f(h)
Поскольку справа линейное отображение, то по определению оно и является дифференциалом [Dx0f]. Мы уже видели примеры таких ситуаций выше, когда рассматривали отображения умножения на матрицу слева или справа. Если f — (скалярная) линейная функция, то она представляется в виде ⟨a,v⟩ для некоторого вектора a — он и будет градиентом f.
Линейность производной. Пусть f(x)=λu(x)+μv(x), где λ,μ — скаляры, а u,v — некоторые отображения, тогда
[Dx0f]=λ[Dx0u]+μ[Dx0v]
Попробуйте доказать сами, прежде чем смотреть доказательство.
И всё бы хорошо, да в первом слагаемом v(x) вместо v(x0). Придётся разложить ещё разок:
[Dx0u](h)⋅v(x)≈
[Dx0u](h)⋅(v(x0)+[Dx0v](h)+o(∣∣h∣∣))=
[Dx0u](h)⋅v(x0)+oˉˉ(∣∣h∣∣)
Это же правило сработает и для скалярного произведения:
[Dx0⟨u,v⟩]=⟨[Dx0u],v⟩+⟨u,[Dx0v]⟩
В этом нетрудно убедиться, повторив доказательство или заметив, что в доказательстве мы пользовались лишь дистрибутивностью (= билинейностью) умножения.
Производная сложной функции. Пусть f(x)=u(v(x)). Тогда
Здесь Dv(x0)u — дифференциал u в точке v(x0), а [Dv(x0)u](…) — это применение отображения [Dv(x0)u] к тому, что в скобках. Итого получаем:
[Dx0u∘v](h)=[Dv(x0)u]([Dx0v](h))
Важный частный случай: дифференцирование перестановочно с линейным отображением. Пусть f(x)=L(v(x)), где L — линейное отображение. Тогда [Dv(x0)L] совпадает с самим L и формула упрощается:
[Dx0L∘v](h)=L([Dx0v](h))
Простые примеры вычисления производной
Вычислим дифференциал и градиент функции f(x)=⟨a,x⟩, где x — вектор-столбец, a — постоянный вектор.
Попробуйте вычислить сами, прежде чем смотреть решение.
Вычислить производную можно непосредственно:
f(x0+h)−f(x0)=⟨a,x0+h⟩−⟨a,x0⟩=⟨a,h⟩
Но можно и воспользоваться формулой дифференциала произведения:
[Dx0⟨a,x⟩](h)=
=⟨[Dx0a](h),x⟩+⟨a,[Dx0x](h)⟩
=⟨0,x⟩+⟨a,h⟩=⟨a,h⟩
Сразу видно, что градиент функции равен a.
Вычислим производную и градиент f(x)=⟨Ax,x⟩, где x — вектор-столбец, A — постоянная матрица.
Попробуйте вычислить сами, прежде чем смотреть решение.
Снова воспользуемся формулой дифференциала произведения:
[Dx0⟨Ax,x⟩](h)=
=⟨[Dx0Ax](h),x0⟩+⟨Ax0,[Dx0x](h)⟩
=⟨Ah,x0⟩+⟨Ax0,h⟩
Чтобы найти градиент, нам надо это выражение представить в виде ⟨?,h⟩. Для этого поменяем местами множители первого произведения и перенесём A в другую сторону (A перенесётся с транспонированием):
⟨ATx0,h⟩+⟨Ax0,h⟩=
=⟨(AT+A)x0,h⟩
Получается, что градиент в точке x0 равен (AT+A)x0.
Вычислим производную обратной матрицы: f(X)=X−1, где X — квадратная матрица.
Попробуйте вычислить сами, прежде чем смотреть решение.
Рассмотрим равенство I=X⋅X−1=I и продифференцируем его:
0=[DX0(X⋅X−1)](H)=
=[DX0X](H)⋅X0−1+X0⋅[DX0X−1](H)
Отсюда уже легко выражается
[DX0X−1](H)=−X0−1⋅[DX0X](H)⋅X0−1
Осталось подставить [DX0X](H)=H, но запомните и предыдущую формулу, она нам пригодится.
Вычислим градиент определителя: f(X)=det(X), где X — квадратная матрица.
Попробуйте вычислить сами, прежде чем смотреть решение.
В предыдущих примерах мы изо всех сил старались не писать матричных элементов, но сейчас, увы, придётся. Градиент функции состоит из её частных производных: ∇x0f=(∂xij∂f)i,j. Попробуем вычислить ∂xij∂f. Для этого разложим определитель по i-й строке:
det(X)=k∑xik⋅(−1)i+kMik,
где Mik — это определитель подматрицы, полученной из исходной выбрасыванием i-й строки и k-го столбца. Теперь мы видим, что определитель линеен по переменной xij, причём коэффициент при ней равен ⋅(−1)i+kMik. Таким образом,
∂xij∂f=(−1)i+kMik
Чтобы записать матрицу, составленную из таких определителей, покороче, вспомним, что
X−1=det(X)1((−1)i+jMji)i,j
Обратите внимание на переставленные индексы i и j (отмечены красным). Но всё равно похоже! Получается, что
∇x0f=det(X)⋅X−T,
где X−T — это более короткая запись для (X−1)T.
Вычислим градиент функции f(x)=∣∣Ax−b∣∣2. С этой функцией мы ещё встретимся, когда будем обсуждать задачу линейной регрессии.
Попробуйте вычислить сами, прежде чем смотреть решение.
Распишем квадрат модуля в виде скалярного произведения:
∣∣Ax−b∣∣2=⟨Ax−b,Ax−b⟩
Применим формулу дифференциала произведения и воспользуемся симметричностью скалярного произведения:
Попробуйте вычислить сами, прежде чем смотреть решение.
Вспомним формулу производной сложной функции:
[DX0u∘v](H)=[Dv(X0)u]([DX0v](H))
и посмотрим, как её тут можно применить. В роли функции u у нас логарифм:
u(y)=log(u),[Dy0u](s)=y10⋅s,
а в роли v — определитель:
v(X)=det(X),[Dy0v](H)=⟨det(X0)⋅X0−T,H⟩,
где под скалярным произведением двух матриц понимается, как обычно,
⟨A,B⟩=i,j∑aijbij=tr(ATB)
Подставим это всё в формулу произведения сложной функции:
[DX0u∘v](H)=det(X)1⋅⟨det(X)⋅X−T,H⟩=
=⟨det(X)1⋅det(X)⋅X−T,H⟩=⟨X0−T,H⟩
Отсюда сразу видим, что
∇X0f=X0−T
Вычислим градиент функции f(X)=tr(AXTX).
Попробуйте вычислить сами, прежде чем смотреть решение.
Воспользуемся тем, что след — это линейное отображение (и значит, перестановочно с дифференцированием), а также правилом дифференцирования сложной функции:
[DX0f](H)=tr([DX0AXTX](H))=
=tr(A⋅[DX0XT](H)⋅X0+AX0T[DX0X](H))=
=tr(AHTX0+AX0TH)
Чтобы найти градиент, мы должны представить это выражение в виде ⟨?,H⟩, что в случае матриц переписывается, как мы уже хорошо знаем, в виде tr(?T⋅H)=tr(?⋅HT). Воспользуемся тем, что под знаком следа можно транспонировать и переставлять множители по циклу:
…=tr(AHTX0)+tr(AX0TH)=
=tr(X0AHT)+tr(HTX0AT)=
=tr(X0AHT)+tr(X0ATHT)=
=tr((X0A+X0AT)HT)
Стало быть,
∇X0f=X0A+X0AT
Вычислим градиент функции f(X)=det(AX−1B).
Подумайте, почему мы не можем расписать определитель в виде произведения определителей
Расписать у нас может не получиться из-за того, что A и B могут быть не квадратными, и тогда у них нет определителей и представить исходный определитель в виде произведения невозможно.
Воспользуемся правилом дифференцирования сложной функции для u(Y)=det(Y), v(X)=AX−1B. А для этого сначала вспомним, какие дифференциалы у них самих. С функцией u всё просто:
[DY0u](S)=⟨det(Y0)Y0−T,S⟩=
=tr(det(Y0)Y0−1S)
Функция v сама является сложной, но, к счастью, множители A и B выносятся из-под знака дифференциала, а дифференцировать обратную матрицу мы уже умеем:
[DX0v](H)=−AX0−1HX0−1B
С учётом этого получаем:
[DX0f](H)=[Dv(X0)u]([DX0v](H))=
=tr(det(AX0−1B)(AX0−1B)−1(−AX0−1HX0−1B))
=tr(−det(AX0−1B)(AX0−1B)−1AX0−1HX0−1B)
Чтобы найти градиент, мы должны, как обычно, представить это выражение в виде tr(?T⋅H).
…=tr(−det(AX0−1B)X0−1B(AX0−1B)−1AX0−1H)
Стало быть,
∇X0f=(−det(AX0−1B)X0−1B(AX0−1B)−1AX0−1)T=
=−det(AX0−1B)X0−TAT(AX0−1B)−TBTX0−T
Вторая производная
Рассмотрим теперь не первые два, а первые три члена ряда Тейлора:
где [Dx02f](h,h) — второй дифференциал, квадратичная форма, в которую мы объединили все члены второй степени.
Вопрос на подумать. Докажите, что второй дифференциал является дифференциалом первого, то есть
[Dx0[Dx0f](h1)](h2)=[Dx02f](h1,h2)
Зависит ли выражение справа от порядка h1 и h2?
Этот факт позволяет вычислять второй дифференциал не с помощью приращений, а повторным дифференцированием производной.
Вторая производная может оказаться полезной при реализации методов второго порядка или же для проверки того, является ли критическая точка (то есть точка, в которой градиент обращается в ноль) точкой минимума или точкой максимума. Напомним, что квадратичная форма q(h) называется положительно определённой (соответственно, отрицательно определённой), если q(h)⩾0 (соответственно, q(h)⩽0) для всех h, причём q(h)=0 только при h=0.
Теорема. Пусть функция f:Rm→R имеет непрерывные частные производные второго порядка ∂xi∂xj∂2f в окрестности точки x0, причём ∇x0f=0. Тогда точка x0 является точкой минимума функции, если квадратичная форма Dx02f положительно определена, и точкой максимума, если она отрицательно определена.
Если мы смогли записать матрицу квадратичной формы второго дифференциала, то мы можем проверить её на положительную или отрицательную определённость с помощью критерия Сильвестра.
Примеры вычисления и использования второй производной
Рассмотрим задачу минимизации f(x)=∣∣Ax−b∣∣2 по переменной x, где A — матрица с линейно независимыми столбцами. Выше мы уже нашли градиент этой функции; он был равен ∇x0f=2AT(Ax−b). Мы можем заподозрить, что минимум достигается в точке, где градиент обращается в ноль: x∗=(ATA)−1ATb. Отметим, что обратная матрица существует, так как rk(ATA)=rkA, а столбцы A по условию линейно независимы и, следовательно, rk(ATA) равен размеру этой матрицы. Но действительно ли эта точка является точкой минимума? Давайте оставим в стороне другие соображения (например, геометрические, о которых мы упомянем в параграфе про линейные модели) и проверим аналитически. Для этого мы должны вычислить второй дифференциал функции f(x)=∣∣Ax−b∣∣2.
Попробуйте вычислить сами, прежде чем смотреть решение.
Вспомним, что
[Dx0∣∣Ax−b∣∣2](h1)=⟨2AT(Ax0−b),h1⟩
Продифференцируем снова. Скалярное произведение — это линейная функция, поэтому можно занести дифференцирование внутрь:
Мы нашли квадратичную форму второго дифференциала; она, оказывается, не зависит от точки (впрочем, логично: исходная функция была второй степени по x, так что вторая производная должна быть константой). Чтобы показать, что x∗ действительно является точкой минимума, достаточно проверить, что эта квадратичная форма положительно определена.
Попробуйте сделать это сами, прежде чем смотреть решение.
Хорошо знакомый с линейной алгеброй читатель сразу скажет, что матрица ATA положительно определена для матрицы A с линейно независимыми столбцами. Но всё же давайте докажем это явно. Имеем hTATAh=(Ah)TAh=∣∣Ah∣∣2⩾0. Это выражение равно нулю тогда и только тогда, когда Ah=0. Последнее является однородной системой уравнений на h, ранг которой равен числу переменных, так что она имеет лишь нулевое решение h=0.
Докажем, что функция f(X)=logdet(X) является выпуклой вверх на множестве симметричных, положительно определённых матриц. Для этого мы должны проверить, что в любой точке квадратичная форма её дифференциала отрицательно определена. Для начала вычислим эту квадратичную форму.
Попробуйте сделать это сами, прежде чем смотреть решение.
Выше мы уже нашли дифференциал этой функции:
[DX0logdet(X)](H1)=⟨X0−T,H1⟩
Продифференцируем снова:
[DX0⟨X−T,H1⟩](H2)=⟨[Dx0X−T](H2),h1⟩=
=⟨−X0−1H2X0−1,H1⟩
Чтобы доказать требуемое в условии, мы должны проверить следующее: что для любой симметричной матрицы X0 и для любого симметричного (чтобы не выйти из пространства симметричных матриц) приращения H=0 имеем
[DX02logdet(X)](H,H)<0
Покажем это явно.
Так как X0 — симметричная, положительно определённая матрица, у неё есть симметричный и положительно определённый квадратный корень: X0=X01/2⋅X01/2=X01/2⋅(X01/2)T. Тогда